library("spectralmc")
# get logprior,
# syntax: lwidth = prior.Gamma(shape, rate) means that lwidth has gamma prior with given shape and rate
logprior <- LogPrior(
amp = prior.Positive(),
lwidth = prior.Gamma(1.01, 0.01), #shape rate
gwidth = prior.Gamma(1.01, 0.01), #shape rate
pos = prior.Unif(0, 1000), # min, max
a = prior.Norm(0, 0.00001), # mean, sd
b = prior.Unif(-1000, 1000)
)
proposal <- Proposal(
amp = prop.Norm(0.005), #sd
lwidth = prop.Norm(0.008),
gwidth = prop.Norm(0.008), # can also use prop.Unif(radius) for uniform proposal
pos = prop.Norm(0.006),
a = prop.Norm(0.00000005),
b = prop.Norm(0.0002) #b = prop.Unif(0.02)
)
fit_kp <- 3 # number of peaks to fit
noise_sigma <- 0.001 # noise_sigma passed to likelihood
signal_model <- voigt.model # which signal to fit (could also be pseudovoigt or GoM)
# get synthetic data
data <- data.synthetic.stormyclouds(seed = 2021, kp = 3, noise = 0.001)
x <- data$x
y <- data$y
posterior <- Posterior(Loglikeli(signal_model, noise_sigma = noise_sigma), logprior)
start_params <- initial.fit.EM(x, y, fit_kp, seed=2021)
|================================================================================| 100% elap: 0s eta: 0s
res <- chains.mhmc(x, y, start_params, posterior, proposal, iter = 20000, print_bar = FALSE)
|================================================================================| 100% elap: 16s eta: 0s
Interactive plot to see fitting process (burn-in of chain): (Not very interesting if EM start is used.)
plt.fit.chain.interactive(x, y, res$samples, signal_model, step_width = 1000)
Double click on legend to zoom to trace.
plt.traceplot(res$samples, "pos")
plt.traceplot(res$samples, "amp")
plt.traceplot(res$samples, "lwidth")
plt.traceplot(res$samples, "gwidth")
plt.traceplot(res$samples, "a")
plt.traceplot(res$samples, "b")
plt.traceplot(res$samples, "b")
# can also do normal coda plots
library(coda)
samples_matrix <- utils.samples_to_matrix(res$samples) #can also filter variables with which=c("pos", "amp")
chain <- mcmc(samples_matrix)
summary(chain)
Iterations = 1:20000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
pos[1] 2.117e+02 1.000e-02 7.074e-05 9.450e-04
pos[2] 5.365e+02 3.752e-02 2.653e-04 1.345e-02
pos[3] 8.666e+02 2.550e-02 1.803e-04 6.722e-03
gwidth[1] 2.780e+00 4.322e-02 3.056e-04 6.808e-03
gwidth[2] 7.517e+00 3.124e-01 2.209e-03 1.627e-01
gwidth[3] 3.975e+00 1.315e-01 9.297e-04 5.933e-02
lwidth[1] 3.740e-01 2.984e-02 2.110e-04 6.793e-03
lwidth[2] 5.862e+00 2.538e-01 1.795e-03 2.774e-01
lwidth[3] 3.445e+00 7.564e-02 5.349e-04 3.708e-02
amp[1] 3.850e+00 1.910e-02 1.351e-04 2.596e-03
amp[2] 4.630e+00 4.187e-02 2.960e-04 1.497e-02
amp[3] 2.363e+00 2.697e-02 1.907e-04 9.631e-03
a 1.911e-06 2.108e-07 1.491e-09 5.145e-08
b 1.938e+01 9.823e-05 6.946e-07 1.267e-05
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
pos[1] 2.117e+02 2.117e+02 2.117e+02 2.117e+02 2.117e+02
pos[2] 5.364e+02 5.365e+02 5.365e+02 5.366e+02 5.366e+02
pos[3] 8.665e+02 8.665e+02 8.666e+02 8.666e+02 8.666e+02
gwidth[1] 2.748e+00 2.766e+00 2.775e+00 2.785e+00 2.818e+00
gwidth[2] 7.136e+00 7.251e+00 7.489e+00 7.660e+00 8.415e+00
gwidth[3] 3.838e+00 3.914e+00 3.954e+00 3.992e+00 4.466e+00
lwidth[1] 3.127e-01 3.631e-01 3.747e-01 3.907e-01 4.190e-01
lwidth[2] 5.468e+00 5.660e+00 5.850e+00 6.093e+00 6.262e+00
lwidth[3] 3.332e+00 3.380e+00 3.428e+00 3.496e+00 3.588e+00
amp[1] 3.829e+00 3.842e+00 3.848e+00 3.855e+00 3.874e+00
amp[2] 4.562e+00 4.602e+00 4.626e+00 4.655e+00 4.750e+00
amp[3] 2.334e+00 2.346e+00 2.359e+00 2.370e+00 2.454e+00
a 1.195e-06 1.827e-06 1.941e-06 2.035e-06 2.206e-06
b 1.938e+01 1.938e+01 1.938e+01 1.938e+01 1.938e+01
plot(chain)
autocorr.plot(chain)